# -*- coding: utf-8 -*-
"""
Creating a normalised model from the original model

@author: obiri
"""
import mpctools as mpc
import numpy as np
import casadi as cs
import time



##16 states, 18 parameters. 

lb = 0.1 # lower bound of normalized model: lb*xs
ub = 1.5 # upper bound of normalized model: ub*xs
width = ub-lb

xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
        6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 
        2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
        5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 
        3.70349577e+01, 4.99220218e+01])

num_est_param = 5
pn=np.zeros((num_est_param))
pn[0]= 6.59e-10/60  #1.0983333333333334e-11
pn[1]= 1.76
pn[2]= 1560.0
pn[3]= 0.06/60      #0.001
pn[4]= 0.075





p_n = np.array([6.59e-10/60, 1.76, 1560.0, 0.06/60, 0.075  ])  ##USED FOR NORMALISATION  IN MAIN FILE
#p_n=np.array([m_mabx, K_damm, rho, mu_dmax, K_gln  ])
p_n =1.1*p_n

num_fixed_param = 13
q=np.zeros((num_fixed_param))
q[0]= 9.6e-03/60  #fix  #0.00016
q[1]= 0.75
q[2]= 28.48 #fix
q[3]= 171.76
q[4]= 2 #fix
q[5]= 0.45
q[6]= 2.0 #fix
q[7]= 2.6e8
q[8]= 8e8
q[9]= 4 #fix
q[10]= 5e5 #fix
q[11]= 1.244 
q[12]= 4e2

q=1.01*q


# np.random.seed(5)
# # std_dev = 0.00009
# # std_dev = 0.00005
# # q = np.random.normal(loc=q, scale=std_dev)    
# # q[7] = np.random.normal(loc=q[7], scale=1e3)    
# # q[8] = np.random.normal(loc=q[8], scale=1e3)    
# # q[10] = np.random.normal(loc=q[10], scale=1e2)    
# # q[12] = np.random.normal(loc=q[12], scale=0.001)    

# q[0] = np.random.normal(loc=q[0], scale=0.000001)  
# q[1] = np.random.normal(loc=q[1], scale=0.00001) 
# q[2] = np.random.normal(loc=q[2], scale=1e-3) 
# q[3] = np.random.normal(loc=q[3], scale=1e-3)   
# q[4] = np.random.normal(loc=q[4], scale=1e-5) 
# q[5] = np.random.normal(loc=q[5], scale=1e-5) 
# q[6] = np.random.normal(loc=q[6], scale=1e-4) 
# q[7] = np.random.normal(loc=q[7], scale=1e4) 
# q[8] = np.random.normal(loc=q[8], scale=1e4) 
# q[9] = np.random.normal(loc=q[9], scale=1e-3) 
# q[10] = np.random.normal(loc=q[10], scale=1e2) 
# q[11] = np.random.normal(loc=q[11], scale=1e-3) 
# q[12] = np.random.normal(loc=q[12], scale=1e-3) 




#plant model
def reactor_tank_uncert(x_aug, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(x_aug.shape[0])
    
    x = x_aug[0:16]
    p = x_aug[16:]
       

    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    p1 =  p[0]*(width*p_n[0])+lb*p_n[0];   p2 =  p[1]*(width*p_n[1])+lb*p_n[1];       p3 =  p[2]*(width*p_n[2])+lb*p_n[2]
    p4 =  p[3]*(width*p_n[3])+lb*p_n[3];   p5 =  p[4]*(width*p_n[4])+lb*p_n[4]


    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    
    
    # K_damm = p[1]#1.76     # ammonia constant for cell death (mM)
    K_dgln = q[0]#9.6e-03/60  # constant for glutamine degradation (min^(-1))
    K_glc =  q[1]#0.75      # Monod constant for glucose (mM)
    # K_gln =  p[4] #0.075     # Monod constant for glutamine (mM)
    KI_amm = q[2]#28.48    # Monod constant for ammonia (mM)
    KI_lac = q[3]#171.76   # Monod constant for lactate (mM)
    m_glc =  4.9e-14/60   ###--ALWAYS FIXED--##   maintenance coefficients of glucose (mmol/cell/min) 
    n =      q[4]#2             # constant represents the increase of specific death as ammonia concentration increases (-)
    Y_ammgln = q[5]#0.45   # yield of ammonia from glutamine (mmol/mmol)
    Y_lacglc = q[6]#2.0    # yield of lactate from glucose (mmol/mmol)
    Y_Xglc =   q[7]#2.6e8    # yield of cells on glucose (cell/mmol)
    Y_Xgln =   q[8]#8e8      # yield of cells on glutamine (cell/mmol)
    alpha1 =   3.4e-13/60   ###--ALWAYS FIXED--##   #  # constants of glutamine maintenance coefficient (mM L/cell/min)
    alpha2 =   q[9]#4        # constants of glutamine maintenance coefficient (mM)
    # mu_max = (0.0016*x[16] - 0.0308)*5/60 # 0.058 * 5  # maximum specific growth rate (min^(-1))
    mu_max = (0.0016*x15 - 0.0308)*5/60    ##--ALWAYS FIXED--##   # 0.058 * 5  ##x[16] changed to the scaled value of x16
    #mu_dmax = (-0.0045*x[16] + 0.1682)/60 # 0.06    # maximum specific death rate (min^(-1))
    # mu_dmax = p[3]# 0.06/60
    # m_mabx = p[0]#6.59e-10/60 # constant production of mAb by viable cells (mg/Cell/min) = Qmab_max
    
    # Add new parameters -- Ben
    Delta_H = q[10]#5e5     # 4.4e4# 1e-8
    # rho = p[2]# 1560.0      # density of mixture is assumed to be density of glucose [g/L]
    cp =      q[11]#1.244        # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    U =       q[12]#4e2
    T_in =  37.0    ###--ALWAYS FIXED--##
    
    ####---these parameters were excluded because they are very easy to know, they willnot change even after a very long time----####
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Diameter of the tank (dm)
    Ac= np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
           
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
            
    # ODEs
    f_lim = (x3 / (K_glc + x3)) * (x4 / (p5 + x4))
    f_inh = (KI_lac/ (KI_lac + x5)) * (KI_amm / (KI_amm + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = p4 / (1 + ( p2 / x6) ** n)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / Y_Xglc + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (alpha2 + x4)
    Q_gln = mu / Y_Xgln + m_gln
    dxdt[3] = (-Q_gln * x1 - K_dgln * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = Y_lacglc * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = Y_ammgln * Q_gln
    dxdt[5] = (Q_amm * x1 + K_dgln * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (p1 * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (Delta_H / (p3 * cp)) * mu * x1 * 2.4908084e-19 + (U / (3.00039503e+03 * p3 * cp)) * (Tc - x15) +w[14]) /(width*xs[14]) 
        
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scaled )



def reactor_tank(x_aug, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(x_aug.shape[0])
    
    x = x_aug[0:16]
    p = x_aug[16:]
       
    #new state indexing
    # Xv1 = x[0]   # concentration of viable cells in reactor(cell/L)
    # Xt1 = x[1]   # total concentration of cells in reactor(cell/L)
    # GLC1 = x[2]  # glucose concentration in reactor(mM)
    # GLN1 = x[3]  # glutamine concentration in reactor(mM)
    # LAC1 = x[4]  # lactate concentration in reactor(mM)
    # AMM1 = x[5]  # ammonia concentration in reactor(mM)
    # mAb1 = x[6]  # mAb concentration in reactor(mg/L)
    # Xv2 = x[7]    # concentration of viable cells in separator(cell/L)
    # Xt2 = x[8]    # total concentration of cells in separator(cell/L)
    # GLC2 = x[9]  # glucose concentration in separator(mM)
    # GLN2 = x[10]  # glutamine concentration in separator(mM)
    # LAC2 = x[11]  # lactate concentration in separator(mM)
    # AMM2 = x[12]  # ammonia concentration in separator(mM)
    # mAb2 = x[13]  # mAb concentration in separator(mg/L)
    # # temperature as a state -- Ben
    # T = x[14]     #temperature of bioreactor mixture  (d C)
    # # buffer tank
    # c = x[15]  # Concentration of mAb in the tank (mg/L)


    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    # p1 =  p[0]*(width*p_n[0])+lb*p_n[0];   p2 =  p[1]*(width*p_n[1])+lb*p_n[1];     p3 =  p[2]*(width*p_n[2])+lb*p_n[2];
    # p4 =  p[3]*(width*p_n[3])+lb*p_n[3];   p5 =  p[4]*(width*p_n[4])+lb*p_n[4];     p6 =  p[5]*(width*p_n[5])+lb*p_n[5];    
    # p7 =  p[6]*(width*p_n[6])+lb*p_n[6];   p8 =  p[7]*(width*p_n[7])+lb*p_n[7];     p9 =  p[8]*(width*p_n[8])+lb*p_n[8]; 
    # p10 = p[9]*(width*p_n[9])+lb*p_n[9];   p11 = p[10]*(width*p_n[10])+lb*p_n[10];  p12 = p[11]*(width*p_n[11])+lb*p_n[11]

    p1 =  p[0]*(width*p_n[0])+lb*p_n[0];   p2 =  p[1]*(width*p_n[1])+lb*p_n[1];       p3 =  p[2]*(width*p_n[2])+lb*p_n[2]
    p4 =  p[3]*(width*p_n[3])+lb*p_n[3];   p5 =  p[4]*(width*p_n[4])+lb*p_n[4]
    

    #inputs remained the same as defined above
    # x8 = xs[7]; 
    # x16 = xs[15];
    # x18 = xs[17];
    
    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    

    #  #parameters - one of the parameters is state dependent so the parameters had to be brought below the normalised variables
    # K_damm = p[1]#1.76     # ammonia constant for cell death (mM)
    K_dgln = 9.6e-03/60  # constant for glutamine degradation (min^(-1))
    K_glc = 0.75      # Monod constant for glucose (mM)
    # K_gln =  p[4] #0.075     # Monod constant for glutamine (mM)
    KI_amm = 28.48    # Monod constant for ammonia (mM)
    KI_lac = 171.76   # Monod constant for lactate (mM)
    m_glc = 4.9e-14/60   # maintenance coefficients of glucose (mmol/cell/min)
    n =   2             # constant represents the increase of specific death as ammonia concentration increases (-)
    Y_ammgln = 0.45   # yield of ammonia from glutamine (mmol/mmol)
    Y_lacglc = 2.0    # yield of lactate from glucose (mmol/mmol)
    Y_Xglc = 2.6e8    # yield of cells on glucose (cell/mmol)
    Y_Xgln = 8e8      # yield of cells on glutamine (cell/mmol)
    alpha1 = 3.4e-13/60  # constants of glutamine maintenance coefficient (mM L/cell/min)
    alpha2 = 4        # constants of glutamine maintenance coefficient (mM)
    # mu_max = (0.0016*x[16] - 0.0308)*5/60 # 0.058 * 5  # maximum specific growth rate (min^(-1))
    mu_max = (0.0016*x15 - 0.0308)*5/60 # 0.058 * 5  ##x[16] changed to the scaled value of x16
    #mu_dmax = (-0.0045*x[16] + 0.1682)/60 # 0.06    # maximum specific death rate (min^(-1))
    # mu_dmax = p[3]# 0.06/60
    # m_mabx = p[0]#6.59e-10/60 # constant production of mAb by viable cells (mg/Cell/min) = Qmab_max
    
    # Add new parameters -- Ben
    Delta_H = 5e5     # 4.4e4# 1e-8
    # rho = p[2]# 1560.0      # density of mixture is assumed to be density of glucose [g/L]
    cp = 1.244        # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    U = 4e2
    T_in = 37.0
    
    ####---these parameters were excluded because they are very easy to know, they willnot change even after a very long time----####
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Diameter of the tank (dm)
    Ac= np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
    
    
    
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
    
    
    
    # ODEs
    f_lim = (x3 / (K_glc + x3)) * (x4 / (p5 + x4))
    f_inh = (KI_lac/ (KI_lac + x5)) * (KI_amm / (KI_amm + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = p4 / (1 + ( p2 / x6) ** n)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / Y_Xglc + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (alpha2 + x4)
    Q_gln = mu / Y_Xgln + m_gln
    dxdt[3] = (-Q_gln * x1 - K_dgln * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = Y_lacglc * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = Y_ammgln * Q_gln
    dxdt[5] = (Q_amm * x1 + K_dgln * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (p1 * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (Delta_H / (p3 * cp)) * mu * x1 * 2.4908084e-19 + (U / (3.00039503e+03 * p3 * cp)) * (Tc - x15) +w[14]) /(width*xs[14]) 
    
    
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scale



#NB:Parameters are not measured and hence are not part of the C matrix
Ny = 7
Nx = 16
Nx_aug = 16 + 5

def measurement_xy_model(x_aug):            
    x = x_aug[0:16]
    
    C_matrix = np.zeros((Ny, Nx))
    C_matrix[0,7]=1
    C_matrix[1,8]=1
    C_matrix[2,9]=1
    C_matrix[3,10]=1
    C_matrix[4,11]=1
    C_matrix[5,12]=1
    C_matrix[6,15]=1   
    y = np.dot(C_matrix, x)
    return y



